Kako bi se sadržak celokupnog genoma nekog organizma uspešno sekvencirao (pročitao), potrebno je DNK podeliti u sitnije fragmente sečenjem pomoću restrikcionih enzima. Dobijeni fragmenata se zatim čitaju odabranom tehnikom sekvenciranja (npr. Illumina ili 454). Dobijene kratke sekvence se slažu odgovarajućim redosledom, čime se dobija sekvenca polaznog genoma. Ovaj postupak sastavljanja sekvenci naziva se asembliranje.
Kratke sekvence dobijene iz sekvencera (eng reads) dele se na kraće, preklapajuće sekvence jednakih dužina, k-mere (eng k-mers), gde $k$ predstavlja dužinu kratke sekvence.
Primer:
sekvenca: ATCGTG
k: 3
k-meri (3-meri):
ATC
TCG
CGT
GTGSledeći korak je konstrukcija grafa. Jedna opcija je konstrukcija grafa čiji su čvorovi k-meri dok grane povezuju preklapajuće k-mere. Rekonstrukcija polazne sekvence ogledala bi se u pronalaženju puta koji prolazi kroz sve čvorove grafa tačno jednom (Hamiltonov put), što je NP težak problem. Druga opcija je konstrukcija grafa čije su grane k-meri spajajući čvorove koji predstavljaju prefiks i sufiks grane. Ovim postupkom se problem svodi na pronalaženje puta koji prolazi preko svake grane tačno jednom (Ojlerov put ili Ojlerov ciklus) što je problem polinomijalne složenosti. Rezultujući graf naziva se De Bruijnov graf.
Za početak, potrebno je pripremiti odgovarajuću klasu Graph koja će čuvati graf, definisan listom susedstva, i pomoćne funkcije koje će biti korišćene za pronalaženje puteva
# Primer liste susedstva
adjecency_list = {
'A': ['B','C'],
'B': ['A'],
'C': []
}
Zatim se strukturi Graph dodaje funkcija koja pronalazi Ojlerov ciklus. Da bi u grafu postojao Ojlerov ciklus, polazeći od bilo kog čvora, čvorovi čiji ulazni i izlazni stepeni su različiti se povezuju dodatnom granom kako bi razlika u stepenima bila uklonjena.
import copy
import random
class Graph:
# Konstruktor grafa, koji za argument dobija listu susedstva
def __init__(self, adjacency_list):
self.adjacency_list = adjacency_list
# Metod vraća listu susednih čvorova zadatog čvora v, uključujući potencijalne duplikate
def get_neighbors(self, v):
return self.adjacency_list[v]
# Metod dodaje novog suseda v čvora u
def add_neighbor(self, u, v):
self.adjacency_list[u].append(v)
# Metod računa izlazni stepen čvora v (broj grana koje izlaze iz čvora v)
def out_degree(self, v):
return len(self.adjacency_list[v])
# Metod računa ulazni stepen čvora v (broj grana koje ulaze u čvor v)
def in_degree(self, v):
degree = 0
for w, neighbors in self.adjacency_list.items():
for u in neighbors:
if u == v:
degree += 1
return degree
# Medod vraća sve grane koje ulaze u čvor v
def inbound_edges(self, v):
edges = []
for u, neighbors in self.adjacency_list.items():
for w in neighbors:
if w == v:
edges.append((u,v))
return edges
# Metod vraća sve izlazne grane iz čvora v ka susedima w u obliku uređenog para (v, w)
def outbound_edges(self, v):
return [(v, w) for w in self.get_neighbors(v)]
# Metod vraća sve grane koje postoje u grafu u obliku uređenih parova (u, v)
def all_edges(self):
edges = []
for w in self.adjacency_list:
edges += self.outbound_edges(w)
return edges
# Metod vraća listu čvorova čiji se ulazni stepen razlikuje od izlaznog stepena
def get_non_balanced(self):
non_balanced = []
for w in self.adjacency_list:
in_deg = self.in_degree(w)
out_deg = self.out_degree(w)
if in_deg != out_deg:
non_balanced.append(w)
return non_balanced
# Metod pronalazi dva nebalansirana čvora i povezuje ih granom tako da više ni jedan
# čvor u grafu nema ulazni stepen različit od izlaznog
def close_to_eulerian(self):
non_balanced = self.get_non_balanced()
if len(non_balanced) != 2:
return None
[u, v] = non_balanced
if self.in_degree(u) > self.out_degree(u):
self.adjacency_list[u].append(v)
else:
self.adjacency_list[v].append(u)
# Metod vraća sve neobiđene grane iz čvora v,
# u odnosu na skup svih neobiđenih grana u grafu
def get_unexplored_edges(self, v, unexplored_edges):
return [edge for edge in unexplored_edges if edge[0] == v];
# Metod proverava da li iz zadatog čvora v ima neobiđenih grana
def has_unexplored_edge(self, v, unexplored_edges):
return len(self.get_unexplored_edges(v, unexplored_edges)) > 0
# Metod proverava da li je graf povezan, u kontekstu DFS obilaska
def is_connected(self):
for v in self.adjacency_list:
visited = set([])
stack = [v]
while len(stack) > 0:
v = stack[-1]
visited.add(v)
has_neighbor = False
for w in self.get_neighbors(v):
if w not in visited:
stack.append(w)
has_neighbor = True
break
if not has_neighbor:
stack.pop()
if len(visited) != len(self.adjacency_list.keys()):
return False
return True
# Metod proverava da li je graf prost,
# da li u grafu postoje čvorovi sa ulaznim stepenom
# većim od 1
def is_simple(self):
for v in self.adjacency_list:
if self.in_degree(v) > 1:
return False
return True
# Metod pronalazi čvor u grafu čiji je ulazni stepen
# veći od 1
def get_non_simple_node(self):
for v in self.adjacency_list:
if self.in_degree(v) > 1:
return v
return None
# Metod formira bypass za čvorove u, v, w
def bypass(self, u, v, w):
# Novi čvor dobija oznaku čvora v uz dodatak
# jedinstvenog stringa u obliku slučajno odabranog broja
x = f'{v}-{random.random()}'
self.adjacency_list[x] = []
self.adjacency_list[u].remove(v)
self.adjacency_list[v].remove(w)
self.add_neighbor(u, x)
self.add_neighbor(x, w)
# Metod proverava da li zadati čvor v ima ulazne i izlazne stepene
# jednake jedinici
def is_1_in_1_out(self, v):
return self.in_degree(v) == 1 and self.out_degree(v) == 1
# Metod pronalazi Ojlerov ciklus u grafu
def eulerian_cycle(self):
# Odabira se proizvoljni polazni čvor
current_node = list(self.adjacency_list.keys())[0]
unexplored_edges = self.all_edges()
# Polazni čvor se dodaje u ciklus
cycle = [current_node]
# Dok postoje neobiđene grane u grafu...
while len(unexplored_edges) > 0:
# Od polaznog čvora se obilaze susedi dok iz čvora postoje neobiđene grane...
while self.has_unexplored_edge(current_node, unexplored_edges):
# Odabira se proizvoljna neobiđena grana iz čvora,
unexplored_from_current = self.get_unexplored_edges(current_node, unexplored_edges)
selected_edge = unexplored_from_current[0]
(u, v) = selected_edge
# odabrana grana se uklanja iz liste neobiđenih grana,
unexplored_edges.remove(selected_edge)
# tekući ciklus se proširuje i pretraga se nastavlja
# od susednog čvora
current_node = v
cycle.append(v)
# Pronađeni ciklus se rotira tako da čvor koji se nalazi u ciklusu
# a iz koga postoje neobiđene grane postane čvor od koga ciklus počinje
for i in range(len(cycle)):
v = cycle[i]
if self.has_unexplored_edge(v, unexplored_edges):
cycle = cycle[i:-1] + cycle[:i + 1]
current_node = v
return cycle
# Metod pronalazi sve Ojlerove cikluse u grafu
def all_eulerian_cycles(self):
all_graphs = [self]
# Dok postoji graf sa cvorom
# koji ima ulazni stepen veći od 1
while True:
found = False
non_simple = None
for g in all_graphs:
if not g.is_simple():
found = True
# Pronaći ne-prosti graf
non_simple = g
break
if not found:
break
# Pronaći čvor sa ulaznim stepenom većim od 1
v = non_simple.get_non_simple_node()
inbound = non_simple.inbound_edges(v)
outbound = non_simple.outbound_edges(v)
# Za svaku kombinaciju ulaznih i izlaznih grana čvora v...
for (u,v) in inbound:
for (v,w) in outbound:
if u == v and v == w:
continue
# Napraviti novi graf sa u, v, w bypass-om
new_graph = Graph(copy.deepcopy(non_simple.adjacency_list))
new_graph.bypass(u, v, w)
# Ako je graf povezan, dodaje se u niz
if new_graph.is_connected():
all_graphs.append(new_graph)
all_graphs.remove(non_simple)
# Izdvajaju se ciklusi u preostalim grafovima
cycles = [g.eulerian_cycle() for g in all_graphs]
unique_cycles = []
# Oznake novonastalih čvorova se isecaju
# do originalnih oznaka čvorova ("v-12345" => "v")
cleaned_cycles = [[v.split('-')[0] for v in cycle] for cycle in cycles]
# Izbacuju se duplikati (izdvajaju jedinstveni ciklusi)
for cycle in cleaned_cycles:
if cycle not in unique_cycles:
unique_cycles.append(cycle)
return unique_cycles
# Metod pronalazi maksimalne, neprekidne putanje u grafu
def maximal_non_branching_paths(self):
paths = []
visited = set([])
# Za svaki čvor v koji nije "1-in-1-out"...
for v in self.adjacency_list:
if not self.is_1_in_1_out(v):
if self.out_degree(v) > 0:
# Od svake izlazne grane iz čvora v...
for (v, w) in self.outbound_edges(v):
non_branching_path = [v, w]
# Beleži se put koji polazi od odabrane grane
# prateći samo "1-in-1-out" čvorove
while self.is_1_in_1_out(w):
visited.add(w)
(w, u) = self.outbound_edges(w)[0]
non_branching_path.append(u)
w = u
paths.append(non_branching_path)
# Pronalaženje izolovanih ciklusa koji nisu obiđeni
# prethodnim obilaskom
for v in self.adjacency_list:
if self.is_1_in_1_out(v) and v not in visited:
visited.add(v)
non_branching_path = [v]
while True:
outbound = self.outbound_edges(v)
if len(outbound) == 0:
break
[(v, w)] = outbound
while w not in visited:
visited.add(w)
non_branching_path.append(w)
outbound = self.outbound_edges(w)
if len(outbound) > 0:
[(w, u)] = outbound
w = u
else:
break
break
return paths
# Primer pronalaženja Ojlerovog ciklusa u proizvoljnom grafu
g = Graph({
'A': ['B'],
'B': ['C', 'E'],
'C': ['D'],
'D': ['A', 'B'],
'E': ['G', 'F'],
'F': ['D','E'],
'G': ['F']
})
g.close_to_eulerian()
g.eulerian_cycle()
Konačno, od sekvence (ili sekvenci) je potrebno konstruisati De Bruijnov graf koji će naslediti sve osobine grafa ali će listu susedstva konstruisati pomoću k-mera.
class DeBruijn(Graph):
def __init__(self, k, reads):
kmers = []
adjacency_list = {}
for read in reads:
n = len(read)
for i in range(0, n - k + 1):
kmer = read[i: i + k]
# u => v
# kmer = ATG
# u = AT, v = TG
u = kmer[:-1]
v = kmer[1:]
if u not in adjacency_list:
adjacency_list[u] = [v]
else:
adjacency_list[u].append(v)
if v not in adjacency_list:
adjacency_list[v] = []
self.adjacency_list = adjacency_list
# Primer konstrukcije De Bruijnovog grafa i pronalaženja Ojlerovog ciklusa
reads = ['TAATGCCATGGGATGTT']
k = 3
dg = DeBruijn(k, reads)
dg.adjacency_list
dg.close_to_eulerian()
cycle = dg.eulerian_cycle()
cycle
# Rekonstrukcija sekvence od pronađenog Ojlerovog ciklusa
# Prvi čvor u ciklusu
seq = cycle[0]
# Od svakog slecećeg čvora uzima se samo poslednji karakter.
# Prolazi se kroz sve čvorove u ciklusu osim poslednjeg, koji zatvara ciklus
for i in range(1,len(cycle) - 1):
seq += cycle[i][-1]
seq
# Pronalaženje svih Ojlerovih ciklusa u DeBruijnovom grafu
dg.all_eulerian_cycles()
dg = DeBruijn(k, reads)
dg.maximal_non_branching_paths()